E 


immmmoE  ^KTEmms  »  ' 

A®i£Dved  te  giUD-iic  r©l®as^ 

,  I>l»sai?3aitos  aaii*2it®4 


1 9960607 


COATING  SELECTION  PROGRAM 


Theory 


DEPARTME?: 
PLASTICS  TECHMIC; 


T  OF  DEFENSE 
,L  £y7iL,ijAT!Of4  C: 


FICiinMNY  ARSEfiAL,  DOVER.  J»- 


by  Frederick  A.  Costello,  Thomas  P.  Harper,  and  Barbara  Aston 


Prepared  by 

GENERAL  ELECTRIC  COMPANY 

Philadelphia,  Pa. 

for 


mic  QUAJOT  mSPEOiaij 


NATIONAL  AERONAUTICS  AND  SPACE  ADMINISTRATION 


WASHINGTON,  D.  C. 


APRIL  196 


TfflS  DOCUMENT  IS  BEST 
QUALITY  AVAILABLE.  THE  COPY 
FURNISHED  TO  DTIC  CONTAINED 
A  SIGNIFICANT  NUMBER  OF 
PAGES  WHICH  DO  NOT 
REPRODUCE  LEGIBLY. 


NASA  CR-1041 


COATING  SELECTION  PROGRAM 


Theory 


By  Frederick  A,  Costello,  Thomas  P.  Harper,  and  Barbara  Aston 


Distribution  of  this  report  is  provided  in  the  interest  of 
information  exchange.  Responsibility  for  the  contents 
resides  in  the  author  or  organization  that  prepared  it. 


Issued  by  Originator  as  Document  No.  65SD526 


Prepared  under  Contract  No.  NASw-960  by 
GENERAL  ELECTRIC  CO.  '  • 
Philadelphia,  Pa. 

for 


NATIONAL  AERONAUTICS  AND  SPACE  ADMINISTRATION 


19960607  090 


RECEIVED  NEW  REPORT,  PAGE  9-9  IS 


BLANK 


JULY  2,  1996 


ABSTRACT 


A  rational  and  direct  method  has  been  developed  for  selecting  the  optical  coating 
pattern  for  the  external ^jjrface  of  a  spacecraft,  such  that  the  spacecraft  tempera- 
tures  are  as  close  as  possible  to  the  midpoint  of  their  preselected  rangea  The 
'tempgrature  control  is  maintained  passively  by  radiation  and  conduction,  using 
no  active  control  devices.  The  complete  rangeTh  mission  envirofiffieffiBs  is  con¬ 
sidered  in  the  optimization  procedure. 

The  selection  technique  has  been  programmed  for  use  on  the  GE  600  Series, 

IBM  7094,  or  any  other  computer  that  uses  the  standard  IBM  Fortran  IV  system. 


FORWORD 


This  document  is  a  three-part  final  report  on  the  Coating  Selection  Program. 

Part  I  describes  the  theory  and  basis  for  selection  of:  (1)  the  iteration  scheme 
used  to  solve  the  heat  balance  equations;  and  (2)  the  optimization  scheme. 

Part  n  describes  the  computer  program,  including  the  details  of  each  subroutine 
and  the  details  of  input  and  output.  Part  n  therefore  includes  the  user's  manual. 
Part  HI  presents  the  Program  Listing.  Parts  I,  II  and  III  have  been  revised 
(Revision  A)  under  a  contract  extension  to  incorporate  descriptions  of  the  first 
two  month’s  usage,  as  well  as  a  new  Program  Listing  with  improvements  as 
found  advisable  during  these  two  months. 

The  work  reported  herein  was  sponsored  by  the  National  Aeronautics  and  Space 
Administration  and  was  monitored  by  Mr.  Conrad  Mook,  of  NASA-Headquarters, 
and  Mr.  Robert  Kidwell,  Jr.,  of  NASA's  Goddard  Space  Flight  Center. 

The  chief  contributors  to  the  work  reported  were  Mr.  Frederick  A.  Costello, 
Engineer,  who  developed  the  techniques,  Mr.  Thomas  P.  Harper,  Analyst,  who 
converted  the  techniques  to  computer  form,  both  of  the  Re-entry  Systems  Depart¬ 
ment's  Thermodynamics  Technology  Component,  and  Miss  Barbara  Aston,  Pro¬ 
grammer-Analyst,  from  the  Spacecraft  Department,  who  programmed  the  tech¬ 
nique  for  computer  usage. 
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1.0  INTRODUCTION 


The  use  of  optical  coatings  to  control  the  temperatures  of  satellites  has  been 
exploited  since  the  first  successful  orbital  flight  in  1958  (Explorer).  Since  then, 
coating  materials  have  been  developed  to  the  point  where  coatings  are  available 
that  give  any  desired  emittance  value  between  0. 1  and  0. 9  for  any  desired  absorptance 
value  between  0. 1  and  0. 9.  What  has  lagged,  however,  is  the  development  of  a  sys¬ 
tematic  approach  to  selecting  coating  patterns  for  the  external  surface  of  the  vehicle, 
such  that  the  satellite  temperatures  are  passively  maintained  as  close  as  possible 
to  the  optimum  temperature.  It  was  in  response  to  this  developmental  need  that  the 
present  work  was  undertaken.  The  work  discussed  in  this  report  constitutes  a 
generalization  of  the  work  performed  by  Costello,  Harper,  and  Cline 

For  the  uninitiated,  an  example  may  prove  useful  in  illustrating  the  effect  of  proper 
coating -pattern  selection.  Figure  1  shows  the  environmental  conditions  for  a  typical 
satellite,  as  well  as  the  optimum  coating  pattern  and  resulting  temperatures.  It  is 
seen  that  using  two  different  coatings  rather  than  one  narrows  the  temperature 
excursions  from  ±7F^  to  ±5F®  from  the  desired  70^F.  This  result  is  not  as 
dramatic  as  those  obtained  for  more  complicated  shapes.  Further  improvement 
is  still  possible,  but  the  point  is  adequately  illustrated. 

The  solution  to  simple  coating  selection  problems,  as  shown  in  Figure  1,  is  indeed 
difficult,  but  the  coniplexity  and  difficulty  increase  rapidly  as  the  number  of  surfaces 
increases  and  as  the  number  of  critical  components  increases.  For  realistic  satel¬ 
lite  designs,  it  is  usually  necessary  to  include  the  entire  conduction  and  radiation 
heat  transfer  networks;  consequently,  the  selection  process  is  obscured  by  the 
various  interactions  and  by  the  number  of  degrees  of  freedom  (twice  the  number  of 
external  surfaces).  Such  complexity  demands  the  use  of  a  high-speed  computational 
device.  From  the  beginning,  then,  the  present  work  has  been  directed  toward  the 
development  of  an  IBM  7094  digital-computer  program  that  would  assist  the  designer 
in  the  coating  selection  process. 

In  the  following  sections,  the  coating-selection  problem  is  formulated  in  mathe¬ 
matically  precise  terms.  The  mathematical  development  of  the  minimization 
scheme  is  presented  in  Sections  3.0,  4.0  and  5.0.  Applications  are  considered 
in  Section  6.0.  The  work  is  summarized  in  Section  7.0  and  future  areas  for 
research  are  cited  in  Section  8.0. 

It  should  be  recognized  that  operating  experience  is  important  to  the  efficient 
usage  of  a  computer  program.  The  present  program  has,  as  all  programs  do,  a 
variety  of  input  constants  that  affect  such  things  as  iteration  procedures  and  opti¬ 
mization  sequences.  An  extention  to  the  original  contract  enabled  a  study  of  the 
input  constants.  The  study  was  conducted  by  analyzing  several  vehicle  designs, 
three  of  which  are  shown  in  section  6.  The  operating  experience  obtained  from 
the  study  along  with  recommended  procedures  for  use  of  the  program  is 
summarized  in  the  User’s  Manual.  It  is  felt  that  the  study,  permitted  by  the  con¬ 
tract  extention,  has  enhanced  significantly  the  benefits  of  the  program. 
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(POSITION  A) 
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POSITION  B) 
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TTTTTj 


Orbit 


Side  1 

Solar  Angle  I  Solar  Flux 


140.0 


Heat  Fluxes 


121.0 


T  =  70°F  is  desired 


Single  Coatini 


(140  +  50. 5)  a  +  (70)e  =  2  e  ct  T^ 
(121  +  43.5)  a  +  (70)  e  =  2  e  cr  T 


Best  Solution 


0.9  0.8  63-77°F 


Double  Coating 


140  ffj  +  50.5  02  +  '2  ""  ^^1  ^2^  ^ 

121  “  j  +  43.  5  o  2  +  ^  2  "  ^^1  ^2^ 

Best  Solution 

“1  “2  ^1  ^  I - ^^-5—1 

0.4  0.5  0.345  0.5  65-75  F 


Hr-Ft‘ 


Coating  Limit 
0.1  ^  a  s  0.9 


0. 1  ^  e  ^0.9 


Figure  1.  Example  of  Two-Coating  Advantages 
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2.0  PROBLEM  STATEMENT 


As  intimated  in  the  Introduction,  the  purpose  of  the  present  work  was  to  develop 
a  computer  program  that  would  determine  the  optical  coating  pattern  that  minimizes 
temperature  excursions  from  some  preselected  temperatures  for  all  possible 
environmental  conditions.  More  precisely,  we  define  the  excursion  parameter. 


th  th 

Where  Ty^  is  the  temperature  of  the  i  network  element  (node)  in  the  k  "  orbit, 

and  Tjy  and  Tjj^  are  the  upper  and  lower  allowable  temperature  limits  of  the  i^^ 

node^  We  wish  to  find  the  surface  coating  properties  of  infra-red  emittance,  €  , 
and  solar  absorptance,  a ,  such  that 

max  5.,  =  min  (2.2) 

i,k 

subject  to  the  constraints  that  a  and  (  are  within  a  preselected  allowable  range. 

From  the  above  definitions,  it  is  obvious  that  (-1/2)  is  the  absolute  minimum  value 
of  and  that  ^  ^  0  implies  that  all  nodes  are  within  their  allowable  temperature 

limits. 

Physically,  the  above  formulation  Implies  that  the  coating  system  is  selected  to 
minimize  the  temperature  excursion  of  the  node  that  is  furthest  from  the  midpoint 
of  its  allowable  temperature  range. 


Other  Possible  Criteria 

The  above  criterion.  Equation  (2.2),  was  selected  from  a  number  of  alternates, 
the  most  interesting  of  which  were 


y 


1 


i,k 


(2.3) 


y 


2 


.a  .  +  0!  .  e.)  =  min 
11  11 


i 


(2.4) 


subject  to  the  additional  constraints 


^0 


Equation  (2.3)  was  the  most  convenient  from  a  mathematical  viewpoint,  but  could 
be  a  minimum  when  one  temperature  was  so  far  out  of  range  as  to  make  the  solution 
ridiculous. 
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Eqviation  (2.4)  resembles  a  linear  programming  formulation  and  is  physically  the 
most  satisfying.  However,  the  equations  are  significantly  non-linear,  and  the 
optimum  may  lie  at  large  distances  from  the  vertices  of  the  restraint  volume. 

Time-Dependent  Solutions 

hiherent  in  the  above  problem  statement  is  the  assumption  that  the  critical  tempera 
ture  can  be  adequately  calculated  using  time-averaged  environmental  heat  fluxes. 
This  is  an  important  and  restrictive  assumption,  complicated  by  the  fact  that  the 
time  constants  of  each  node  are  different  and  vary  non-lin early  with  the  node 
temperature.  In  Section  8.0,  further  study  is  recommended  in  this  area.  In  the 
meantime,  however,  the  average-flux  procedure  should  serve  as  a  useful  guide 
in  selecting  the  optimum  coating  pattern.  The  average-flux  analysis  can  be  made 
quite  accurate  if  the  emittance  is  kept  low  and/or  the  node  heat  capacities  kept 
high,  so  that  all  the  time  constants  are  significantly  greater  than,  say,  the  orbit 
period. 

Development  of  Solution 

The  solution  to  Equation  (2.2)  can  be  conveniently,  but  not  completely,  divided 
into  two  steps:  (1)  solving  for  the  T., 's,  given  a  set  of  a's  and  e ’s;  and  (2)  select¬ 
ing  new  values  of  a  and  f  to  reduce  S .  The  first  problem  is  discussed  in 
Section  3.0;  the  second,  in  Sections  4.0  and  5.0. 
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3.  0  SOLVING  THE  HEAT- BALANCE  EQUATIONS 


The  node  temperatures  are  governed  by  the  heat-balance  equations: 

conduction  radiation 

interchange  interchange 


SK..  (T.,  -  T., )  + 

1]  Ik  jr 


EffR..(T,  -T.,)  = 
1]  Ik  ]k^ 


=  “i  %i  ^  (^Sli  ^ik  -  "i  "  Tik>  S: 


incident  incident 
solar  albedo 

flux  flux 


incident  re-emission  internal 
earth  heat 

flux  generation 


absorbed  solar  energy  net  absorbed  IR  energy 

where 

i,  j  =  node  numbers  =  1,  2,  . , . ,  N  +  S 
k  =  orbit,  or  time-interval,  number  =  1,  2,  . . . ,  q 
and  where  the  nomenclature  is  described  in  Appendix  D  (Section  9. 0)< 
Equation  (3. 1)  can  be  written  in  the  more  compact  form 

-  4  4 

EA  T.  +  '51^  M..  T.  =  a  .  S.  +  e  .  E.  -  a .  O'  T.  +  Q. 

13  3  2-f  ^  ^  ^  1  1  1  1 


i  =  1,2,  ...,  N 

where  for  convenience  we  have  dropped  the  k  (orbit  number)  subscript.  In 
discussing  solutions  to  Equation  (3.2),  it  is  frequently  necessary  to  examine  the 
symmetrically  linearized  form: 

yL  T-c.  (3.3) 


where 


-  K..  -  CT  R..  (T.  +  T.)  (T^.  +  T.^)  .  /  . 

1]  1]  ^  1  '  1  ]  '  1  7^  3 


L..  = 
13 


-  V  L..  +  4  a.  e  .  CT  T: 

11  1  1  ”  i 
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;;  4, 


C.  =  Ci  (Kgi  s. K^.A.) (K^  .  E. 3  a.  a  )  Kq  Q. 


where  once  again  k  has  been  dropped. 


Many  methods  for  solving  the  linear  system,  Equation  (3. 3),  have  been  devised. 

A  convenient  and  up-to-date  summary  is  given  by  Fox (3).  These  methods  have 
been  adapted  to  the  present  set  of  equations  and  the  speeds  of  convergence  com¬ 
pared.  In  addition,  two  new  matrix-inversion  methods  have  been  examined.  In 
the  following  paragraphs,  each  method  is  described.  All  methods  make  use  of  the 
residual  vector,  r^^: 


•i  ^ 


(3.5) 


When  the  exact  solution  is  obtained,  r^  =  0  for  all  values  of  "i". 
Gauss-Seidel  Methods 

Applied  to  the  linear  system.  Equation  (3. 3),  the  extrapolated  Gauss-Seidel 
method  gives; 

C,-E 

rj,  '  ^  i<i  ^  ^ _ 3 

i  -  L.. 


T  (n+1)  ^  T  . 

i  oi  i 


(n) 


a 


n 


n 


In  terms  of  the  residual  vector,  this  may  be  written 

r.(>^)x  V  T  (Symmetrically 

^3  Linearized 


T. 

1 


(n+1)  ^  T  + 
i 


ill 


o^n^ii 


Gauss- 

Seidel) 


(3.6) 


The  factor  is  the  extrapolation  (  “^  ^  1)  or  interpolation  (  >  1)  factor. 

Equation  (3.6)  was  used  directly  for  the  present  non-linear  case,  with  the 
residual  vector  defined  by  Equation  (3.  5). 

In  the  non- stationary  iteration  process,  is  governed  by  the  equation 

^(n+1)  _  ^(n) 

“n  ~  ^  ^(n)  _  rp(n-l) 
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However,  from  the  e3q)erience  of  Kaplan'  %  it  is  necessary  to  restrict  a  between 
1/2  and  1  to  obtain  convergence  for  the  heat- balance  equations.  For  the  Resent 
case,  a  was  taken  as  1. 0,  a  constant. 

Some  saving  in  computation  can  be  realized  if  the  linearization  process  applied 
to  Equation(3. 2)  takes  a  Taylor-series  form: 

^ .  r  3  ~ 


Ek..  T.  +y'  M..  -  y  (A  .. 

1]  ]  Z-/  ij  ]  ['  1] 


+  4M..  T.  )  T.  -  3  M..T 
1]  1  J  1]  ] 


In  this  case,  the  iteration  formula  becomes 

-  S  [a..  +  4M..  T.^  1  1 

^(n+1)  ^(n)  ,  ^  i<i  L  ^  ^  ^  JL  3  3  J 

^  ^  [  ^ii  ^  ^ii  ^i^  ] 


Maximum  Rate  of  Descent 


Applied  to  Equation  (3. 3),  the  maximum  rate  of  descent  method  described  by  Fox 

r  E  r/”)  r/-) 

T/n^l)  =  T/n)  +  (3.8) 


ri  13  1  ] 


This  was  carried  over  directly,  using  Equation  (3.  5)  to  calculate  r^. 
Coniugate  Gradient  Method 


Fox  gives  for  the  linear  system 

FE 


^(n+1)^  (n)  _2 

i  i  V 


(11)  (n) 
w.  '  r. 

]  1 

(n)  (n] 

L..  w.  '  w.  ^ 
1]  1  3 


where 


n  =  0 


n  s  1 


Once  again,  using  Equation  (3.  5)  to  define  r..  Equation  (3. 9)  can  be  used  in  the 
non-linear  case.  3 
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Newton-Raphson  Method 


The  Newton-Raphson  method  uses  Equation  (3.  5)  directly,  taking 


Setting  r.  equal  to  zero  gives 


j"^  (n+1)  _  ^  (n)j 


(n+1)  _  „  (n)  V  N  r 

i  i  Y  ^  /X 

3  /  \  (n) 

where  ■)  N. .  >  is  the  inverse  of  the  matrix  (  -tit  1  where 
\  Ijj  \5Tj/ 


ar.  \  (n) 


A..  -  4  M..  T  . 
1]  1]  ] 


3  - 

4  e.  a.  CT  T  -  A  -  4  M..  T.  i  =  j 
111  1]  ij  ] 


This  method  requires  a  matrix  inversion  at  each  iteration  and  was  found  slower 
than  the  matrix- inversion  method  described  next. 

One  way  of  alleviating  the  difficulty  of  invertir^  a  matrix  each  time  is  to  use  an 
approximate  inverse  in  the  form 


(E  +  e)  =  [e  (  1  +  E"^  c)] 

=  [l  -  E"^  ej"^  E“^ 

=  (I  -  E"^  f  )  E"^ 

Then  the  expression  for  T^^  ^  becomes,  for  E^  =  ( ^  N  =  E 


l  rrt  _  rr  (®)  TvT  - 


(n+l)  _  ip  (n)  _  jq-  J.  W  ^ 

L  i  ij  j 


-  N  r  N  f  N  r^"^  (Expanded 

ij  j  D  ^ij  jk  k  1  1  Inverse) 

j  j,k,l 


where 


[(T.^"^)^  -  (T.^®V 

1]  L'  ]  '  J  '  J 


a.  a  +  M.  J  •  [(T.^^'b^- (T.^V]  i  =  j 
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Matrix  Inversion  Methods 


A  direct  Matrix  Inversion  method  of  solving  Equation  (3. 2)  can  be  obtained  simply 
by  writii^; 


i 


X.  +  a  .  s.  +  f.  E.  -  a.  cr 
1  ]  ]  ]  ]  3  3  3  3 


(3. 10) 


where  differs  from  Ay  in  that  a  convergence  factor,  X.  ,  has  been  added, 
according  to  the  convergence  criteria  of  Appendix  B  (Section^  9. 0).  Note  that 
Equation  (3. 10)  converges  quite  rapidly  vhen  the  radiation  terms  are  small. 

Using  the  residual  vector.  Equation  (3. 10)  can  be  written  in  the  form 

(Direct  Matrix  Inversion)  (3. 11) 


where 


J'..>  =  (A’.. 
Ui  )  13 


-1 


which  is  evaluated  once  per  solution. 

One  improvement  to  Equation  (3. 11)  can  be  obtained  by  writing  Equation  (3. 2) 
in  the  form 


A..  +  M.. 
13  13 


=  A..  +  M..  i  i 

13  13  ^ 


A"..  = 
13 


V  A"..  +  4  e.^o)  ^  ^  (T.^°y  +  X.  s  A..  +  M.. 
/-«  1]  1  1  1  '  1  11  11 

3 


1  =  3 


and 


L(A..  T.  +  M..  T.'^  +  M  ..  =  a,  S.  +  e,  E,  -  €  ct  a  T  ^  + 

'  1]  3  13  3  13  3  ^  1  11  1  11 


Then 


^U+l)  =  2  j-  «j  Sj  +  "  3  Ej  -  V  T,  ^ik  ^ 

j  1]  [  k 

-  L 

k 


In  terms  of  the  residual  vector 


(n+1)  _  .j.  (n)  ^  ^  jtt  r  (•^)  (Linearized  Matrix  Inversion)  (3. 12) 

1  -  i  j  «  j 

J"  is  evaluated  once  per  solution. 

i] 

Comparison  of  Methods 

It  is  beyond  the  scope  of  this  development  program  to  develop  ^e  system  of 
Loremsrequiredto  compare  the  rates  of  convergence  of  the  above  me*^  of 

solviiK  the  heat-balance  equations.  Such  a  task  is  rendered  difficult 
heterogeneous  character  of  the  iteration  schemes.  However,  it  is  practical 
importont  to  evaluate  the  schemes  numerically  for  laical  satellite-type  hea  - 
balance  equations.  The  methods  are  compared  m  Tables  I  and  H. 

There  are  several  parameters  that  must  be  defined  before  the  significance  of  the 
SlowiS  numericaf  comparisons  can  be  understood.  Convergence  is  assumed  to 

be  completed  if 

I T  ^  DT  i  =  l,  2,  •••9  N 


j  -  r.^”^  s  DR 


i  =  1,  2,  . . . ,  N 


In  addition,  ttose  methods  that  use  T,  the  linearizing  temperature,  update  T  (i.  e. , 
let  T  =  every  time 

I  T  -  T.^”^  s  DTR  i  =  1,  2,  . . . ,  N 

I  i  ^ 

Those  methods  using  the  convergence  parameter,  X,  calculate  the  value  of  X 

Tj^aX  '^SOL- 
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TABLE  L  SOLUTION  TIMES  FOR  ITERATION  SCHEMES 


Iteration  Scheme 

A 

PROBLEM 

B  C 

1. 

Direct  Matrix  Inversion 

2.2 

1. 0  1, 77 

2. 

Linearized  ]V&trix  Inversion 

1.8 

1. 1  1. 27 

3. 

Expanded  Inverse 

4.2 

2. 0  3.  32 

4. 

Linearized  Gauss-Seidel 

1.3 

1. 3  1. 79 

5. 

Symmetric  Gauss-Seidel 

1.8 

1.7  2.63 

6. 

MRD 

2.7 

2.7  3.19 

7. 

Conjugate  Gradient 

1.6 

1.2  1.42 

NOTES 

1. 

Problems  Used  Were: 

A.  Pure  Radiation  (h^A  ^  1, 0  BTU/hr°F, 

between  adjacent  nose) 

B.  Pure  conduction,  except  external  surfaces  which  used  pure  radiation 

(hj,A-  1.0) 

C.  Combined  conductiop/ radiation;  h^A 

kA/x=  1. 

0 

In  all  cases,  8,  16,  and  64  node  problems  were  considered.  The 
geometry  was  a  cube  broken  into  equally  sized  cubes.  The  8  node  case 
used  a  2  X  2  X  2  configuration;  the  16,  a  2  x  2  x  4;  and  the  64,  a  4  x  4  x  4. 

2. 

In  general,  DTR  =  DR  =  DT  =  1. 0  for  the  above  cases.  A  20%  to  30% 
improvement  was  made  in  all  methods  but  (1)  and  (2)  by  setting  DTR  =  0. 0; 
that  is,  by  not  correcting  the  original  linearization. 

3. 

Methods  1,  2,  and  3  require  specification  of  Tiji^^vto  assure  convergence. 

The  correlated  time  is  good  for  Tj^j^  ^  1,  3  Tsql*  ^  20%  increase  in 

ruiming  time  was  ejqjerienced  for  T,,. ^  =  1. 5  T„, . 

MAX  SOL 

4. 

See  Table  n  for  actual  nmning  times  and  additional  notes® 
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TABLE  n.  TIMING  DATA  FOR  ITERATION  SCHEMES 


Problem  Nodes  At 


Iteration  Scheme  (Times  In 


Millisec 


(Nodes)"* "  (AT)''* 

5  6  7 


2.6  1.68  1.40  3.08  1.13 


2.4  1.61 
2.4  2.34 
15.6  1.45 
15.6  1.55 
15.6  2.25 

2.4  1.12 
2. 4  . 563 


64  15. 6  . 446 

8  2. 4  1.  68 

16  2. 4  1.  69 

64  2. 4  2. 19 

8  5. 9  1. 56 

16  5.  9  1.  50 

64  5.  9  2. 17 

8  15.  6  1. 58 

16  15.6  1.44 

64  15.  6  2. 13 

8  15.6  .526 

16  15.6  .378 

64  15. 6  . 375 


.28  3.37  1.21 

.  58  4. 77  1. 63 

.  18  2. 89  0. 92 

.  13  3.06  0.95 

.63  4.56  1.45 

.84  1.96  1.12 

.723  1.45  1.12 
.633  1.49  1.61 
.79  1. 58  0. 92 

.  68  1. 36  1. 02 

.  66  1.  5  1. 36 

L.  40  3. 08  1.  68 

L.29  2.98  1.77 

L.  41  4. 24  2. 42 

L.  17  2. 92  1.  56 

L.  17  2.90  1.56 

L.40  4.39  3.24 

L.  18  2. 76  1.  32 

L.  10  2.87  1.51 

1.  39  3. 96  2. 05 

.789  1.32  4.48] 
.642  1.28  6.65 
.673  1.39  7.63' 


2.24 

2.81 

4.13 

2. 14 
2.40 
3.  36 
1. 72 
2.04 
2. 84 

6.97] 

12.3 

1.91' 


8.  6® 

15.2^^^ 


NOTES; 


1.  For  these  runs,  DT  =  1. 0  and  DR  =  0. 5,  but  the  converged  solutions 
were  3  to  7°F  in  error. 

2.  For  these  runs,  DT  =  1. 0  and  DR  =  0. 5,  but  the  runs  were  terminated  for 
being  unsuccessful  after  100  iterations.  Solutions  are  normally  obtained 
in  10-20  iterations. 

3.  For  this  run,  DT  =  DR  =  1. 0,  but  the  converged  solution  was  about  80°F 
in  error. 

4.  AT  is  the  change  (  in  %)  in  equilibrium  temperature,  based  on  the  initial 
temperature  in  *^R,  from  the  initial  temperature  to  the  exact  final  temper¬ 
ature  (TgQ^). 

5.  The  geometry  used  was  a  cube  broken  into  equally  sized  cubes. 
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Examination  of  Tables  I  and  n  reveals  that  Method  (2)  is  the  fastest  in  nearly  every 
case.  One  notable  exception  is  the  conduction  dominated  case  (B),  where  Method 
(1)  was  faster.  The  reasonfor  this  is  that  (2)  used  too  conservative  a  stability  factor 
for  these  particular  cases.  If  the  normal  values  had  been  used  (see  Appendix  B), 
Method  (2)  would  have  shown  the  approximately  20%  improvement  over  Method  (1) 
that  is  evident  in  the  other  cases  analyzed. 

One  other  case  was  studied  that  is  not  reported  on  the  tables.  This  was  a  67-node 
model  of  an  ablative  heat  shield  for  a  satellite/ reentry  vehicle.  For  this  problem, 
due  to  the  peculiarities  of  the  conduction  matrix,  Ay,  only  Methods  (1)  and  (2)  con¬ 
verged,  the  times  correlating  well  with  those  shown'^on  the  tables.  The  other  methods 
were  no  more  than  1%  of  the  way  toward  a  solution  after  100  iterations,  indicating 
that  their  convei^ence  times  would  be  a  factor  of  one  thousand  greater  than  the  times 
indicated  on  Tables  I  and  H. 

The  choice  of  methods,  therefore,  is  clearly  in  favor  of  Method  (2). 
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4.0  GEOMETRY  OF  THE  g  SURFACES 

Of  paramount  importance  in  the  development  of  any  optimization  scheme  is  fore¬ 
knowledge  of  the  geometry  of  the  criterion  surface.  For  example,  if  there  is 
more  than  one  minimum  point,  a  method  must  be  devised  for  first  finding  one 
minimum,  then  search  for  the  second,  and  so  on,  to  determine  which  point  is  the 
absolute  minimum.  In  the  following  paragraphs,  unimodality  and  other  properties 
of  the  -surface  are  examined. 

Monotonic 

The  first  important  characteristic  of  /3  is  that  it  is  composed  of  segments  of  mono¬ 
tonic  surfaces.  That  it  is  composed  of  segments  of  surfaces  can  be  seen  from  the 
definition  of  ^  (see  Equation  2. 2).  That  is  composed  of  two  monotonic  branches 
can  be  seen  as  follows: 

The  heat-balance  equation  can  be  written: 

E  AyTj  +E  t  K^jOj 

3  3 

(4.1) 

4 

This  equation  can  be  linearized  in  T  according  to  the  scheme 


k-1 


wiiere 


i  ^  j 

I  ^li  ■'■  Z  (^ii  ■'■  /  "^iP 

\  i^j  ^T.**  Mi 

(4.3) 

i  =  j 

Note  that  since  N.,.  ^  0  for  i  ^  j  and  N„  =  -S  N.., 
1]  11  1] 

{Aj}”'  =“ » 

(4.4) 

-1  .... 

This  property  of  |N.j|  holds  regardless  of  the  value  of  the  T/s,  since  Tj^  ^  0. 

Equation  (4. 2)  carries  on  the  right-hand  side  one  term  involving  T^.  This  term 

vanishes  if  (a)  there  is  no  conduction,  or  (b)  there  are  no  constant-temperature 
nodes  to  which  heat  is  conducted. 


Equation  (4. 2)  gives 


N  a,p4 


A  ST. 

e.aa.  ^ 


r  1 


(4.  5) 


3T? 


N  hrf 


N 


(4.6) 


N4S 


(4.7) 

N+S  4 

Em..t: 

1]  1 

j=Nfl 


Since |Ny  +  has  the  same  properties  of  {n^^, 


Therefore,  from  Equation  (4,  6),  T  is  monotonically  increasing  with  a. .  To 
determine  how  T7  varies  with  define  ^ 


Then 


Y.  =  cra.T.  -  K^.E. 
1  11  ^11 


3Y.  3Tr 

=  oa. 


^/N..  \  3Y. 


.  I  oa.  1  1, 


)N  N  N-fS 

Y,  -  U,(K,,S.  .  .  g  ^  E,  .  -  I  Z  -  i,) 


j=N^l 


(4.10) 


yT.  MRHS). 


Therefore,  from  Equation  (4.  9) 


— ^  =  P  Y 


'C  1 —  X 


(4. 11a) 


Since  |p„|  ^  0,  Equation  (4. 11)  indicates  that  Yj^  increases  or  decreases  with 
dependii^  on  the  signs  of  the  right-hand- sides. 

E  (KgjS.  +  ^  0  for  all  values  of  "i",  then  Equation  (4. 6)  shows  that  no  T^  - 

surface  has  a  horizontal  point.  Ji  (%%  +  KaiAi)  =  0  for  all  "i".  Equation  (4. 11a) 
must  be  examined,  along  with  the  defihition  of  (RHS)j,  Equation  (4. 10). 


—  =  ^ 

,  1=1  L  3-1  J 


N+S 

I  E  5.(T.  -  T  )  -  E  M  .tJ 
^j=N+l  ^  ^  i=Nfl  ^ 

(4.  lib) 


Within  the  accuracy  of  the  linearization,  the  expression  for  (RHS)j  is  a  constant, 
independent  of  the  values  of  a  and  e  .  Its  sign,  therefore,  will  not  change  over  a 

given  Yk  (or  Tk^)  surface.  Since  {Pk^}  is  positive,  also  will  not  change 

sign  on  the  Yk  surface. 

4 

The  conclusion  to  be  drawn  from  the  above  is  that  the  Yk>  und  therefore  Tk 
and  Tk  (see  Equation  (4.  8)),  have  no  extremes  inside  the  boundaries  of  any  given 
Yjj  surface. 

A  similar  analysis  linearizing  T^  into  T  leads  to  the  same  conclusion.  Trying 
to  solve  the  non-linear  case  directly  is  frustrated  by  the  lack  of  an  explicit 
expression  for  Yk  in  Equation  (4. 9). 

Consider  next  a  composite  surface,  p,  defined  by 

p  =  max  p 
i,k  ^ 


where 


p.,=  max 


=  max 


=  max 


4  4 

"Vd-  "4 

4  4 


[(Cj  -  Cj)  ;  (Cj  -  C, 


Define 


Cl  Yik  -  C2 


^ikL"C2 


Cl^ik 


Now  the  Piui  and  /3iku  have  the  same  properties  of  the  Yk's  analyzed  above. 
EquationsT^  6)  and  (4. 11)  imply  that  the  minimum  value  of  ^must  occur  at  the 
intersection  of  at  least  two  jS  surfaces. 


Unimodality  of  jS 


The  next  objective  is  to  show  that  the  preceding  properties  lead  to  the  conclusion 
that  there  can  be  only  one  relative  minimum  point.  Suppose  there  were  two,  say 
i3  j  and  ^  2*  Then  around  each  of  these  points 

6p>0  (4.12) 


Therefore  there  is  a  relative  maximum,  jS  o,  lying  along  a  continuous  but  broken 
path  following  the  axes  between  Points  1  and  2.  Then  around 

6  P  ^  0  (4. 13) 


Consider  a  line  through  Point  3  parallel  to  the  x.  axis.  A  variation  in  jS  along 
the  axis  is  designated  fl  Now  the  condition  (4.  i3)  implies  that 


t  -  MSTANCE  ALCWG 
PATH  1  -  2 


Although  /Sis  discontinuous  at  Point  3,  the  surfaces  that  intersect  there  are  sections 
of  infinite,  continuous,  and  monotonic  surfaces,  so  that 


Therefore,  if  6/3.^  >  0,  <  0,  contradicting  Equation(4. 13).  Therefore  there 

cannot  be  two  discrete  minima.  However,  if  6  =  0,  there  may  be  an  infinite 

number  of  minima  connected  by  a  continuous  path. 


Concavity  of  g 

A  surface  is  simple  or  concave  if  for  any  two  points  on  the  surface,  p  j,  and 

Xj  +  (1  -  X)x2j  +(  1  -  x)  jsj^XgJ 

K  we  consider  a  line  along  the  c  .  axis,  Equation  (4. 11)  shows  that  p is  either 
convex  or  concave.  Therefore  j^has  segments  that  are  concave  and  raier  segments 
that  are  convex,  pis  therefore  not  strictly  concave  or  convex. 

Figure  2  shows  a  typical  p  surface,  from  which  the  jg surface  is  constructed. 
Figure  3  shows  a  typical  map  for  a  one-node/one-external-surface  problem. 
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Figure  2.  Section  of  a  Typical  8^^  Surface 


^-6 


Figure  3.  Map  of  |3  Surface  for  One-Node  Case 


5.  0  OPTIMIZATION  SCHEMES 


Several  optimization  schemes  were  considered  for  use  in  the  present  proeram 
These  included: 

(a)  Linear  Programming 

(b)  Pattern  Search 

(c)  Miguel's  Poor-Man's  Ricfee  Follower 

(d)  Variational  Methods 

(e)  Hill- climbing  (Maximum  Rate  of  Descent,  or  MRD) 

Methods  (b)  and  (c),  which  are  described  in  Wilde's  book:  "Optimum  Seeking 
Methods, '  can  fail  to  find  the  minimum  due  to  the  extremely  sharp  ridges  (dis¬ 
continuous  derivatives)  characteristic  of  the  -  surfaces. 

Method  (d),  used  in  many  trajectory  optimization  schemes,  is  frustrated  by  the 
fact  that  has  discontinuous  derivatives.  It  is  possible,  however,  to  redefine 
/3  and  consider  the  independent  variables  to  be  a.,  and  T.,  subject  to  the 
constraints  ^  ^  ^ 


^  T 

<  T 

iL  i 

“  iU 

<  a 

S  0( 

iL  i 

iU 

^  ^ 

£  € 

iL  i 

iU 

This  leads  to  a  rapid  solution  if  the  solution  satisfies  these  limits  within  the 
constraints  of  the  heat  balance  equation. 


The  most  severe  disadvantage  of  the  variational  scheme  is  the  number  of  compu¬ 
tations  required  if  the  constraints  cannot  be  met.  For  example,  a  three-node 
problem,  each  with  independent  values  of  «  ,  e  ,  and  T,  has  N  =  9  independent 
variables.  Considerirg  both  upper  and  lower  bounds  on  the  variables,  the  heat- 
balance  equation  must  be  solved  S  times,  where 


S  =  2(N-1) 


=  13180 

In  the  MRD  scheme,  S  is  on  the  order  of  100  to  200.  Since  most  of  the  computation 
time  is  spent  solving  the  heat-balance  equation,  the  variational  schenie  need  not  be 
considered  further. 
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Method  (a),  although  quite  powerful  in  linear  problems,  depends  on  the  accuracy  o 
linearization  for  nonlinear  problems.  It  is  anticipated,  on  the  basis  of  past  ex¬ 
perience,  that  some  of  the  external- surface  temperatures  will  range  from  say 
300°R  to  500^^  due  to  orbit-to-orbit  changes  in  the  relative  position  of  the  sun. 
Linearization,  therefore  appears  dangerous.  In  addition,  if  there  no  feas^le 
solution,  a  rather  extensive  sensitivity  analysis  would  be  necessary  to  determme 
what  must  b6  don©  to  obtain  a  feasible  solution. 

Suppose,  however,  that  these  objections  were  overcome.  The  linear-programming 
time  per  iteration  would  be  on  the  order  of 


Time  (  nsec)  =  1000  p  +  61  (1  +  p)  +  400  p  q  r 

(2  rn  +  2o)  (2  rp  +  4q)^  x  10^ 

+  -  I2 


where 

N  =  number  of  nodes 
jp  =  number  of  orbits 
q  =  number  of  external  nodes 
r  =  number  of  critical  nodes 

For  reasonably  sized  problems^  the  last  term  dominates.  In  fact,  one  finds  ihe 
time  nearly  proportional  to  (rp)^.  The  following  table  gives  some  time  estimates 

for  q  equal  to  10. 


p 

Time  (sec) 

A- 

5 

25 

10 

130 

5 

200 

10 

1060 

5 

1570 

10 

8450 

5 

5300 

10 

28500 

For  problems  of  this  general  size,  MRD  would  take  an  estimated  3,600  seconds. 

From  the  foregoing  discussion,  it  may  be  concluded  that  tte  MRD  approach 
most  promising  due  to  the  h^hly  non-linear  character  of  the  equations,  the  long 
runni^  times  for  LP  systems  with  no  feasible  solutions,  and  ^e 
nii^  ttoes  between  LP  and  MRD  for  systems  of  common  complexity  and  with 
feasible  solutions.  Hence,  MRD  was  selected  for  the  present  program. 

The  MRD  method  is  described  in  many  texts^®^.  It  involves  selecting  a  starting 
point  (a  set  of  ai’ s  and  f^'s)  and  calculating  the  criterion  function  {p)  as  well  as 
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the  values  of  (a/S/s  a^)  and  (a  ^/a  at  this  point.  Using  the  derivatives,  the  di¬ 
rection  in  which^-|^  is  a  maximum  can  be  determined  where 


is  a  maximum  can  be  determined  where 


(ds)2  = 


(dep2^(d«.)2 


Then  for  a  given  (as),  A  jS  will  be  a  maximum,  and  j3  will  decrease  more  than  if 
any  other  direction  had  been  used  (provided  As  is  small  enough).  Thus  a  step-by- 
step  procedure  is  used  to  decrease  j3,  where  each  step  is  taken  in  the  local  maximum- 
rate-of-descent  direction. 


MRD  converges  to  the  optimum  very  slowly  when  it  enco\mters  a  sharp  ri<^e,  such 
as  are  characteristic  of  the  p  -surface  defined  by  Equation  (2. 1).  Several  ridge- 
followit^  techniques  are  described  in  the  literature,  but  these  are  frustrated  by  the 
extreme  sharpness  of  the  present  ri^es.  None  use  the  ridges  to  assist  in  finding 
the  solution.  Therefore,  a  method  called  TREND  was  devised  for  the  present  problem. 

In  TREND,  use  is  made  of  the  fact  that  the  jS  ji,  surfaces  are  monotonic,  so  that 
when  a  ridge  is  crossed,  say  from  Point  1  toToint  2,  one  of  the  subscripts  of 
changes  (or  a  switch  from  upper  to  lower  bound  occurs).  If  the  ridge  is  immediately 
re- crossed  to,  say  Point  3,  this  indicates  that  the  points  being  studied  lie  in  a 
valley.  A  vector  extrapolation  from  Point  1  through  Point  3  gives  an  improved 
value  of  /3  (see  Appendix  A).  Figure  4  illustrates  the  process.  The  vector  PI,  P3 
is  in  the  downward  direction  of  the  valley  and  indicates  the  trend  of  the  valley.  In 
many  test  cases,  it  has  been  found  that  TREND  results  in  a  time  saving  of  1.  5;  1  to 
10:1  over  the  conventional  MRD  methods. 


The  logic  diagram  for  the  computer  program  developed  in  the  present  work,  using 
the  TREND  procedures,  is  shown  in  Appendix  C  (Section  9. 0). 


c 


Figure  4.  Illustriition  of  TREND  Method 


6.0  APPLICATIONS 


FICTITIOUS  SATELLITE 


As  an  example  of  the  use  of  the  foregoing  theory,  the  problem  of  optimizing  the 
coating  pattern  of  a  cylindrical,  horizontally  stabilized  earth  satellite  will  be  con¬ 
sidered.  Figures  5  and  6  show  the  details  of  the  satellite  geometry,  the  orbits 
considered,  and  the  arrangement  of  the  coating  patches. 

Table  HI  shows  the  heat  fluxes  (BTU/hr-ft^)  on  each  patch  for  the  range  in  solar 
angles  shown  on  Figure  6.  These  are  the  average  heat  fluxes  over  one  orbit 
period. 

The  details  of  the  input/output  procedures  and  the  dctual  computer  program  are 
shown  in  Section  5  of  Part  n  of  this  final  report.  The  results,  however,  are 
summarized  on  Table  IV.  The  importance  of  the  optimization  process  is  obvious 
on  examination  of  Table  IV.  Even  for  an  optimized  single- coating  pattern,  the 
temperature  excursions  from  the  desired  mean  are  50%  greater  than  can  be  attained 
with  the  multi-patch  coating. 

NASA  EPE-D 

A  second  example  was  supplied  by  Mr.  Robert  Kidwell  of  NASA-GSFC.  This  was 
a  mathematical  model  of  the  NASA  EPE-D  (Explorer  XXVI)  Satellite.  The  nodal 
breakdown  is  shown  on  Figure  7  and  a  photograph  of  the  actual  satellite  is  shown 
on  Figure  8.  The  orbit  of  this  satellite  is  such  that  the  heat  flux  received  from 
the  earth  is  negligible.  The  significant  parameter  is  the  angle  between  the  solar 
ray  and  the  spin  axis.  In  this  case,  instead  of  orbit  number,  the  different  sets  of 
heat  fluxes  represent  different  angles  between  solar  ray  and  spin  axis. 

The  engineers  at  GSFC  had  had  some  experience  with  the  temperature  response 
of  satellites  of  this  configuration.  These  investigators  undertook  a  trial- and- error 
procedure  using  considerable  engineering  judgement  to  select  a  suitable  coating 
pattern  for  the  EPE-D.  Mr.  Kidwell  estimates  that  at  least  4  weeks  and  4  hours 
of  computer  time  were  required  to  achieve  the  desired  results.  Table  V  shows  the 
temperature  ranges  attained,  along  with  the  allowable  ranges.  Only  the  critical 
nodes  are  shown  on  this  table. 

The  challenge  was  to  determine  if  the  computer  program  could  be  used  to  improve 
on  this  design  and  the  design  procedure.  The  Coating  Selection  Program  did  de¬ 
vise  two  better  coating  patterns.  The  results  of  the  one- week  study  are  also 
summarized  on  Table  V.  A  total  of  approximately  40  minutes  of  computer  time 
was  used.  This  shows  that  the  program  can  reduce  costs  by  a  factor  of  four. 

NASA  IMP-C 


A  third  example  was  again  supplied  by  Mr.  Robert  Kidwell  of  NASA-GSFC.  This 
was  a  mathematical  model  of  the  NASA  IMP-C  satellite.  The  nodal  breakdown  is 
shown  in  Figure  9,  and  a  photograph  of  the  actual  satellite  is  shown  in  Figure  10. 
The  Coating  Selection  Program  devised  a  better  coating  pattern  than  the  one  de¬ 
vised  by  the  engineers  at  GSFC.  The  study  using  the  Coating  Selection  Program 
took  one  week  and  50  minutes  of  computer  time  which  also  represents  a  significant 
gain  over  the  original  design  procedure  used  at  GSFC.  The  results  are  summarized 
in  Table  VI. 
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10’ 


Figure  5.  Vehicle  Geometry 


RANGE  IN  SOLAR 
FLUX  ANGLES  =  ±  60° 
FROM  OVERHEAD 


VIEW  AA  -  8  PATCH 
Figure  6.  Orbit  and  Patch  Geometry 
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TABLE  m.  ORBITAL  HEAT  FLUXES 


00 


TABLE  IV.  RESULTS  OF  COATING-OPTIMIZATION  STUDY 


Specifications:  Allowable  «  -  e  Range:  0. 1  s  e  £  0. 9,  0. 1  ^  a  ^  0. 9 

Allowable  Temperature  Range:  ^ 

Skin  Temperatures:  200  R  to  800  R 
Internal  Temperature:  525°R  to  5450R 

Results: 
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TABLE  V:  SUMMARY  OF  EPE-D  DESIGNS 


NODE 

Temperature  Ranges 

Allowable  Limits 

Predicted  Limits 

Using  NASA 
Coatings 

Using  Computer  Solution 

#1 

#2 

#3 

1 

420 

564 

468  - 

537 

438  - 

489 

478  - 

542 

411  - 

492 

6 

456 

- 

600 

471  " 

600 

482  - 

578 

477  - 

572 

432  - 

618 

7 

456 

- 

600 

468  - 

599 

478  - 

583 

475  - 

574 

434  - 

624 

8 

456 

- 

600 

470  - 

578 

479  - 

566 

482  - 

566 

443  - 

611 

17 

492 

- 

546 

495  - 

535 

496  - 

539 

499  - 

540 

486  - 

538 

18 

492 

- 

546 

509  - 

524 

510  - 

529 

501  - 

540 

493  - 

524 

20 

474 

- 

582 

488  - 

564 

497  - 
_ _ 

574 

517  - 

572 

456  - 

580 

NODE 

Coatings  Used  _ 

Allowable  Limits 

By  NASA 

Using  Computer  Solution 

— ii 

#2 

#3 

O' 

e 

a 

e 

a 

c 

a 

e 

1 

.278 

.260 

.350 

.849 

.152 

.118 

.97 

.86 

2 

.170 

.112 

.121 

.032 

.149 

.115 

.391 

.826 

3 

.170 

.112 

.121 

.034 

.153 

.116 

.424 

.81 

4 

.170 

.112 

.123 

.038 

.155 

.118 

.474 

.776 

5 

.320 

.260 

.129 

.040 

.150 

.119 

.510 

.768 

6 

.310 

.240 

.236 

.232 

.167 

.196 

.561 

.457 

7 

CO 

o 

CO 

C<| 

CO 

00 

lO 

00 

.310 

.240 

.254 

.20 

.168 

.199 

.556 

.560 

8 

• 

o 

II 

• 

o 

II 

• 

o 

l] 

• 

o 

11 

.310 

.240 

.256 

.179 

.168 

.198 

.679 

.574 

9 

U/ 

Vl/ 

w 

w 

.970 

.850 

.946 

.853 

,207 

.166 

.902 

.793 

10 

csi 

tH 

CO 

CO 

t- 

05 

in 

CO 

.970 

.850 

.944 

.859 

.146 

.123 

.365 

.667 

11 

o 

s 

o 

& 

o 

II 

a 

d 

II 

a 

,970 

.850 

.941 

.859 

.170 

.147 

.75 

.644 

12 

tH 

CO 

.160 

.110 

.143 

.108 

.163 

.183 

.705 

.599 

13 

o 

*rH 

"S 

.310 

.240 

.123 

.039 

.153 

.123 

.431 

.792 

14 

s 

s 

Si 

.350 

.850 

.123 

.039 

.144 

.116 

.359 

.764 

NOTES: 

1.  Computer  solutions  differed  only  in  their  starting  points,  which  were  as  follows: 


1  0=0. 35,  e=0. 85  for  node  1;  a=0. 97,  e=0. 86  for  nodes  9, 10, 11;  a=0. 12,  e-0. 03 
for  all  other  nodes. 

2  e=ctrO.  12  for  all  nodes 

3  0=0. 97,  €=0. 86  for  node  1;  o=0. 35,  €=0. 85  for  all  other  nodes 

2.  The  third  computer  case  terminated  without  attaining  an  acceptable  solution.  It 
appears  that  if  the  step  size  had  not  been  reduced  so  rapidly,  a  better  solution 
would  have  been  attained. 
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TABLE  VI:  SUMMARY  OF  IMP-C  DESIGN 


)E 


Temperature  Ranges 

Allowable  Limits 

Predicted  Limits 

Using  NASA 

Using  Computer  Solution 

Coatings 

501 

555 

499 

539 

504 

545 

501 

555 

483 

532 

505 

529 

492 

555 

474 

526 

496 

539 

501 

555 

492 

532 

505 

534 

474 

600 

528 

591 

515 

592 

492 

555 

503 

530 

505 

530 

501 

555 

503 

553 

504 

535 

474 

564 

537 

578 

538 

560 

501 

555 

496 

537 

504 

537 

Allowable  Limits 

Coatings  Used 

By  NASA 

Using  Computer  Solution 

ry 

ry 

€ 

.17 

.03 

.613 

.48 

.25 

.85 

.12 

.03 

1.0 

.84 

.628 

.496 

03 

04 

CD 

00 

m 

00 

.18 

.188 

.486 

.387 

«  o 

• 

• 

.35 

.84 

.299 

.704 

II  II 

II 

II 

.21 

.03 

.504 

.364 

w  w 

W 

Vt/ 

.22 

.188 

.247 

.786 

.7 

.81 

.591 

.670 

.256 

.24 

.181 

.415 

.74 

.82 

.724 

.598 

cq  csj 

tH  • 

m 

CM 

.279 

.636 

.244 

.815 

t- 

Oi 

.22 

,188 

.149 

.218 

o 

« 

• 

.87 

.82 

.745 

.621 

II  !l 

II 

II 

.256 

.157 

.133 

.111 

e 

.32 

.188 

.188 

.459 

.35 

.84 

.435 

.718 

7.0  CONCLUSIONS 


The  theoretical  foundation  and  a  detailed  logic  diagram  have  been  developed  from 
which  a  digital  compiiter  program  has  been  devised  that  will  select  a  space-vehicle 
external-coating  pattern  in  such  a  way  that  the  internal  temperatures  of  the  vehicle 
can  be  maintained  as  constant  as  possible  in  a  pvtrely  passive  manner. 

The  program  has  been  demonstrated  on  several  vehicle  designs,  three  of  which  are 
shown  in  Section  6. 0  of  this  report.  It  is  seen  that  the  gain  of  usii^  a  multiple- 
patch- external  coating  amounts  to  a  factor  of  2  in  temperature  variations  from 
orbit  to  orbit.  The  coating  pattern  is  derived  in  a  direct  method  in  a  few  minutes 
time  on  a  digital  computer  with  approximately  2 -microsecond  access  time  (such  as 
an  IBM  7094).  The  added  ejqiense  of  the  enhanced  design  is  seen  to  be  quite  small. 

The  details  of  the  program  will  be  foimd  in  Part  H:  User’s  Manual  and  Program 
Description. 
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8. 0  RECOMMENDATIONS  FOR  FURTHER  WORK 


Time-Dependent  Solutions 

The  present  work  assumes  that  the  critical  temperatures  can  be  calculated  using 
average  heat  fluxes,  that  is,  that  the  component  temperature  deviates  little  from 
the  average  temperature.  This  assumption  is  good  if  the  node  heat  capacitance  is 
high  compared  to  the  equivalent  thermal  conductance  between  the  node  and  its  en¬ 
vironment  so  that  the  time  constant 

(see  Equation  3. 4)  is  much  greater  than  the  time  interval,  t,  for  which  the  average 
heat  fluxes  are  calculated.  If  is  much  greater  than  one  orbit  period,  no  diffi¬ 
culties  arise.  * 

For  thin  external  surfaces,  t  .  is  usually  quite  small,  so  that  temperature  devia- 

tions  from  the  mean  are  large.  It  is  therefore  important  to  determine  the  true 
average  temperatures  and  the  magnitude  of  the  deviation  from  the  mean. 

For  transient  heat  flow.  Equation  (3. 2)  becomes 


where 


A  A 

("Vi  Vi  -Dj(e)  -  yjOTj 


D.(0)  -  f  jEj  +  Qj 


Defining  X  = 


J  X  de,  i 


integration  of  Equation  (8. 2)  gives 


1  o  e 


("Vi  [ Vi 

L  —  J  3 

If  the  temperature  history  is  periodic  of  period  Tp,  then,  if  Sj  =  +  Tp,  the 

heat-capacitance  term  is  zero.  If,  in  addition,  A. .  h  0,  the  average  fluxes  will 

4  —  ” 

give  the  true  average  T  's.  If  A.,  y  0,  then,  since 

*J 

T^ 
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neither  the  true  average  nor  T  are  calculated.  This  is  a  problem  ttat  shouM 
be  considered  in  future  work.  One  approach  would  be  to  devise  form  factors,  y , 
such  that 

=  y 

However,  the  development  of  this  concept  is  beyond  the  scope  of  the  present  work. 
Convergence  of  Iterative  Solution 

The  heat-balance  equation  was  solved  for  several  numerical  examples  usu^  seven 
iterative  schemes.  The  fastest  of  the  methods  was  then  selected  for  the  CSP.  The 
rates  of  convergence  should  be  investigated  and  compared  on  a  more  mathematical 
and  rigorous  basis. 

Sub-optimizations  of  each  iteration  scheme  should  also  be  conducted.  For  example, 
one  could  determine  if  X  could  be  made  smaller  for  iteration  schemes  (1)  aM  (2) 
(Table  I),  or  one  mi^ht  add  X  to  the  L  matrix  of  method  (7)  to  determme  if  this  aids 
convergence  in  such  cases  as  Problem  D  (Table  n). 


"Typical"  Orbits 

In  selecting  the  orbits  used  to  obtain  the  optimum  coating  pattern,  the  usual  approa^ 
is  to  take  the  extremes  in,  say,  ;3-angles  (i.  e. ,  the  angle  brtween  solar  ray  and 
plane),  and  then  take  a  "few"  orbits  in  between.  This  would  give  representative  orbits 
for  auniformly  coated  vehicle.  However,  when  the  coating  pattern  is  optimized  for 
the  particular  orbits  given  in  the  input,  it  becomes  questionable  whether  the  particular 
orbits  still  represent  extreme  cases. 


For  most  spacecraft,  the  environmental  heat  fluxes  can  be  expressed  as  a  function 
of  one  or  two  variables  (/5-angle,  solar  view-angle,  orbit  position,  etc.),  so  a 
solution  to  the  problem  of  properly  selecting  the  orbits  should  be  attainable. 


Program  Enlargement 

The  program  is  now  limited  to  34  nodes.  For  many  satellites  this  is  not  adequate. 
The  program  can  be  enlarged  if  a  UNK  (chain)  procedure  is  used,  so  that  tape  and 
drum  storage  is  used. 


Tolerances 

The  present  program  uses  the  nominal  values  of  a  and  e.  The  effects  of  tolerances 
in  these  values  on  the  space  vehicle  temperatures  must  be  calculated  separately. 
This  extra  calculation  could  be  incorporated  in  the  optimization  program. 


Heaters  and  Shutters 

The  present  program  considers  only  those  vehicles  in  which  no  auxili^y  tempera¬ 
ture  control  device  is  used.  As  such,  it  determines  the  best  control  that  can  be 
attained  without  such  devices.  The  next  step  in  determining  the  optimum  tempera¬ 
ture  control  system  would  be  to  include  heaters  and  shutters  (sometimes  called 
louvres,  vanes,  or  variable-surface-property  device^. 
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9.0  APPENDICES 


Appendix  A,  Mathematical  Analysis  of  the  TREND  Step 

The  mathematical  basis  for  TREND  can  be  stated  as  follows: 

Given  two  intersecting  hypersurfaces,  A  and  B.  Given  also  a  point  ’a'  on  A  and 
a  point  ’b’  on  B.  ’b’  lies  on  the  projection  of  the  maximum-rate-of-descent  line 
of  A  through  ’a*,  the  projection  being  taken  on  the  hyperplane  A  (x.)  =  0,  The 
projected  distance  from  *a’  to  ’b*  is  'As*.  Given  also  a  point  'c*  oK  A  having 
the  same  relationship  to  ’b'  as  ’b’  does  to  'a'. 

We  wish  to  prove  that  A  (c)  ^  A(a).  If  this  is  the  case,  then  B  (d)  ^  B(b),  where 
’d’  is  found  using  the  maximum-rate-of-descent  line  from  c.  This  process  can  be 
continued  imtil  a  minimum  value  of  A  or  B  is  reached.  This  will  prove  that  the 
maximiun-rate-of -descent  method  will  converge  to  the  minimum.  For  TREND, 
we  must  show  also  that  A  (e)  s  A(c),  where  i-c  =  [(c-a)/lc-al  ]  As. 

Consider  first  the  maximum  rate  of  descent.  The  point  b  is  given  by 


bj=a,- 


A  .  (As) 

XI 


where  the  independent  variables  are  taken  as  x..  Also 


Ci  =  b. 


B  .  (As) 

XI  '  ' 


Therefore 


Now 


c.  -  a. 
1  1 


A  .  B  . 

XI  ^  XI 


As 


A(c)  =  A(a)  +  ^  A^(a)  *  (c.  -  a^  + 


=  A(5)  - 


A".  A  .B  . 

XI  ^  XI  XI 


As 
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Therefore 


from  ’a’,  say 


Similarly 


Ui  •  Ug  =  cose 

where  e  is  the  angle  between  Uj  and  Ug.  Therefore 
tJj  *  tfj  +  •  Ug  =  1  +  cose  s  0 


or 


Therefore  A(c)  ^  A(a),  provided  As  is  small  enough  that  the  first  term  of  toe 
Taylor's  expansion  is  representative  of  toe  function  A  (no  more  than  say  1/2%  in 
error).  If  toe  radius  for  which  the  single  term  expansion  is  adequate,  is  designated 
R,  it  is  required  that  |c  -  aj  ^  R,  or 
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2 


(-Tfe'Tl? 


(As)^  ^  (A-1) 


The  maximum  value  of  the  term  inside  the  brackets  is  4,  so  it  is  sufficient  if 


As  ^  J  =  0. 5R 


(A-2) 


In  TREND,  we  wish  to  state  also  that  A(e)  is  A(c),  where 


(c  -  a ) 

e.  -  c.  =  — . . .  (As) 


Now  if  (as)  is  chosen  small  enough 

A(e)  =  A(a)  +  2A^.(a)  (e.  -  a.)  + 


=  A(5)+  1  + 


E(c, -aj2  i 


2^  A  .(a)  (c,  -  a.) 


XI'  '  '  1  1' 


Therefore 


(e.  -  c. 

1 

c.  -  a. 

1  1 


A(i)  -  A(5)  =  — --g-—  ^  X) 


[a(c)  -  A(a) j  ^  0 


which  proves  the  validity  of  TREND,  The  (As)  must  be  chosen  such  that 


e  -  a  £  R 
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or 


E 


(«i  -  a/  “ 


or 


1  + 


A  s 


E  (Cj  -  \ 


E(Ci-/^R=' 


but 


S  (c.  -  a.)^  ^  4  (A  s)^  From  Equation  (A-1) 


1  1 


••  E  (Cj  -  a/  +  2(A  s)  VE(Ci-a/  +  (A  s)^  ^  4(A  sf  +  4(A  s)^  +  (A  s)^  =  9  (As)‘ 


Therefore,  As  must  be  such  that 
As  s:  ^  R  =  0. 333R 


(A-3) 


It  is  seen  by  comparison  between  Equations  (A-2)  and  (A-3)  that  TREND  will  be 
slower  than  MRD  if  the  A  .’s  equal  the  B  /s;  that  is,  if  the  two  steps  of  MRD  are 

approximately  in  the  same  direction.  However,  if  Cj^=  so  that  E(c.  -  a.)  =  e  , 

MRD  will  move  with  net  step  sizes  on  A  of  e  and  convergence  will  be  slow.  TREND, 
under  the  same  circumstances  will  proceed  with  a  net  step  size  of  (e  +  As), 
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Appendix  B.  Convergence  of  the  Inverse-Matrix  Iteration  Methods 

To  analyze  the  convergence  of  the  matrix  inversion  iteration  schemes  (Methods 
1  and  2  of  Table  I),  it  is  necessary  to  Unearize  the  radiation  terms.  Let 


-CT  R..  (T.  +  T.)(Tf  +  T?)  i5^j 
N+S_ 


(B-1) 


Then  Equation  (3. 2)  becomes 


N+S 


where  C.  is  defined  as  in  Equation  (3.4). 
Method  (1)  solves  this  equation  in  the  form 

=  (1  +  X)"^  (C  -  (R  -  X)  T^"^) 


(B-2) 


Convergence  of  the  iteration  process  depends,  then,  on  the  spectral  radius, P  (Z),  of 

rr  _  /A  .  /-.j-  «\ 


Z  s  (A  +  X)  ^  (X  -  R) 
being  less  than  one. 

Note  that 


(B-3) 


X)„  »  -E  A,J 

and 

A.. 

1] 

^  0 

N 

Ry  »  -  E  Ry 
1=1  •* 

and 

«ii 

^  0 

i  j 

p: 

IV 

o 

and 

=  0 

(B-4) 


The  inatrix  (A  +_X)  is  positive  by  virtue  of  Equations  (B-4).  K,  in  addition 
Xy  2:  R„,  (  X  -  R)  is  positive  so  that,  since 


(A+R)  =  (A  +  X)  -  (X  -R)  with  (A  +  R)~^  2;  q. 


(B-5) 


[(A  +  X)  -  (X  -  R)]  is  a  regular  splitting  of  the  matrix  (A  +  R).  Then  by  Theorem 
2.2  of  VargaW 


p(Z)  =  p[(A  +  R)'^X-R)]  <  J 


P  [(A+R)'^  (X-R)] 


(B-6) 
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so  that  convergence  is  assured  if  Xii  =  Rii.  By  Theorem  2. 7  of  Varga,  it  is 
evident  that  Xii  =  Rii  gives  the  minimum  value  of  P  (Z)  and,  therefore,  the  most 
rapid  convergence. 

Relating  the  foregoing  to  the  non-linear  equation  shows  that  a  sufficient  condition 
for  convergence  is  that 


N  +  S  9  9  3 


(B-7) 


and  that  most  rapid  convergence  will  occur  if 


o  9  9  3 

Xji  =  -  E  Rij  (Ti  +  T.)  (Tf  +  T.)  +  4  a.ej<TT. 

iJl 


(B-8) 


where  Tk  and  X^  are  evaluated  at  each  iteration.  This  latter  criteria  can  be 
approximated  by  using  Ti  as  Ti,  solution. 

It  is  evident  that  if  Equation  (B-8)  is  satisfied.  Equation  (B-7)  will,  in  general,  be 
violated.  Although  it  has  not  been  shown  in  a  mathematically  rigorous  f^luon  ex¬ 
perience  with  many  numerical  cases  indicates  that  the  criteria  e3q)ressed  in  Equa¬ 
tion  (B-7)  can  be  relaxed  somewhat  so  that  convergence  will  occur  if 


X..  =  max 


[I  ^^ii  ■  V  " 


(B-9) 


where  R*  =  R-  evaluated  with  T.  and  T.  at  their  maximum  anticipated  values. 
1]  1]  I  J 

Combining  Equations  (B-8)  and  (B-9)  gives  the  following  criteria  on  X.^: 


Xii  =  max  j^J(R*.-A..);E|A..|-A.j;E(-Rij  +  Ay)-Aj.;Ryj  (B-10) 

(Method  1) 

where  R°i  =  R^  evaluated  with  T.  and  Tj  at  their  anticipated  solution  values. 

The  analysis  of  Method  (2)  follows  the  saine  pattern  as  for  Method  (1).  The  counter 
part  of  Equation  (B-10)  is  found  by  replacing  R  in  Equation  (B-10)  by  R  -  R  . 

Xii  -  max  (R*  -  R^  -  A^);ElA^I  -  Aii;E(-Rij  +  Ry  +  Ai*)  -  Aii;  0 

•'  (B-11) 

(Method  2) 
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Appendix  D. 

Symbol 

A 

A 

C 

D 


E 


K 


Qi 


L 


M 


N 


P 


Q 

R 


(RHS) 

S 

T 

(WCp), 

X 


Nomenclature  (Excepted  and  Noted) 

Description 
incident  albedo  flux 
conduction  matrix  term 
defined  by  Equation  (3. 4) 
defined  by  Equation  (8. 2) 
incident  earth  flux 
(arj/3T.)  at  T.  = 
inverse  matrix  of  A 

multiplier  of  incident  albedo  flux  for  node  *i’ 

multiplier  of  incident  earth  flux  for  node  ’i’ 

multiplier  of  internal  heat  generation 

conductance  from  node  'i'  to  node  'j* 

multiplier  of  incident  solar  flux  for  node  ’i* 

defined  by  Equation  (3. 4) 

radiation -matrix  term 

defined  by  Equation  (4. 3) 

inverse  of  N..  +  e.  6.. 

1]  1  1] 

internal  heat  generation 

defined  by  Equation  (9. 2. 1) 

product  of  radiant  interchange  factor  between 
nodes  'i'  ajid  ’i'  and  area  of  node  'i* 

defined  by  Equation  (4, 10) 

incident  solar  flux 

absolute  temperature 

heat  capacitance  of  node  ’i' 

stability  parameter  in  matrix  inversion  methods 


Suggested  Units 
BTU/hr  ft^ 
BTUAr  °F 
BTU/hr 
BTU/hr 
BTU/hr  ft^ 
BTU/hr  °F 
hr  °F/BTU 
ft2 
ft^ 

BTU/watt  hr 
BTU/hr  °F 
ft^ 

BTU/hr  °F 
BTU/hr  ft^  (°R)^ 
BTU/hr  (°R)^ 
hr  (°R)^/BTU 
watts 

BTU/hr  °F 
ft2 

BTU/hr 
BTU/hr  ft^ 

°R 

BTU/°F 
BTUAr  °F 
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V>^ 


Symbol 


Description 


Suggested  Units 


y 

C 

d 

r 

ct 

a 

n 

P 

y 


€ 


0 

X 

a 

r  . 

Cl 

Subscripts 


re -radiating  area  for  node  'i*  (external  area) 
weighting  factor,  Equation  (2. 3) 
weighting  factor,  Equation  (2. 4) 
residual  vector  (Equation  3. 5) 
solar  absorptance 

iteration  scheme  extrapolation  parameter 
temperature  deviation  parameter 
criterion  function.  Equations  (2, 3)  and  (2, 4) 
form  factor,  ^)/(T)^ 

Kronecker  delta  (=1  if  i=i;  =o  if  i?^j) 

infrared  emittance 

time 

interpolation  parameter 

-8 

Stephen-Boltzmann  constant  (0. 1713  x  10  ) 

(WCp)/Li, 


i  usually  node  number 

iL  corresponding  to  lower  limit  on  T^, 

iU  corresponding  to  upper  limit  on 

k  usually  orbit  number 


Superscript 


BTU/hr 


BTU/hr 


hrs 


BTU/hr  ft^  (°R)^ 
hrs 


n 


iteration  number 
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